import numpy as np
import Mixing_rules as MR
from Pure_VLE_TB import T_BOIL_CAL
from EOS_collection import CUBIC_EOS
from scipy.optimize import minimize

R = 8.314

#输入对应态参数
Tc = np.array([[364.8, 408.1]])
Pc = np.array([[4.610e6, 3.648e6]])
w = np.array([[0.148, 0.176]])
Z = np.array([[0.275, 0.281]])
M = np.array([[44, 44]])
Tb = np.array([[322.02, 373.42]])

#输入要计算的状态参数
x1 = np.array([[0.3]])
x2 = 1 - x1
x = np.concatenate((x1, x2))
ynew = x.flatten() * Pc.flatten() / ( np.sum ( x.flatten() * Pc.flatten() ) )
T0 = np.ones_like(x) * 300
T = np.sum(x.flatten() * Tb.flatten())
P = 2026.5e3
n = 10 #迭代循环次数

EOS = CUBIC_EOS(Pc, Tc, Z, w)


def BUBBLE(Tb):

    for i in range(n):
        if i < 2:
            T = Tb[i]

        # 初始化
        ynew = x

        # 计算状态方程参数
        a = EOS.PR_a(T)
        b = EOS.PR_b()

        # 计算液相逸度系数
        aL = MR.a_mix(a, x)
        bL = MR.b_mix(b, x)
        VL = 2 * bL
        VLnew = EOS.PR_V(P, VL, T, aL, bL)
        while np.sum ( np.abs(VLnew - VL) ) / np.sum ( VL ) > 1e-5:
            VL = VLnew
            VLnew = EOS.PR_V(P, VL, T, aL, bL)
        ZML = P * VL / R / T
        xa = MR.xya_mix(a, x)
        phiL = EOS.PR_PHI_MIX(P, VL, T, aL, bL, xa, b, ZML)

        # 求气相的组成
        # 初始化
        tol = 1
        while tol > 1e-5:
            # 计算气相逸度系数
            y = ynew
            aV = MR.a_mix(a, y)
            bV = MR.b_mix(b, y)
            VV = R * T / P
            VVnew = EOS.PR_V(P, VV, T, aV, bV)
            while np.sum ( np.abs(VVnew - VV) ) / np.sum ( VV ) > 1e-5:
                VV = VVnew
                VVnew = EOS.PR_V(P, VV, T, aV, bV)
            # VV=VL，说明此时此温度低于泡点，气相y为0
            if np.abs(VV - VL) < 1e-5 :
                y = np.zeros_like(y)
                break
            ZMV = P * VV / R / T
            ya = MR.xya_mix(a, y)
            phiV = EOS.PR_PHI_MIX(P, VV, T, aV, bV, ya, b, ZMV)
            ynew = x.flatten() * phiL / phiV
            ynew = ynew.reshape(-1, 1)
            #此处，如果选取的温度离泡点较近，容易出现y的值振荡而不收敛，视情况这里取平均值帮助收敛
            ynew = (y + ynew) / 2
            tol = np.abs(np.sum(y - ynew))

        if i == 0:
            F0 = np.sum(y) - 1
            T0 = T
        else:
            F = np.sum(y) - 1
            #割线法迭代
            T1 = T - F * ( T - T0 ) / ( F - F0 )
            T0 = T
            T = T1
            F0 = F
        
    return T

T = BUBBLE(Tb.flatten())
print(f"泡点温度是{T}K")